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INTRODUCTION 


The  rapid  discharge  of  propellant  gas  from  a  weapon  produces  a  strong 
shock  wave  which  propagates  into  the  environment.  Several  shocks,  contact 
surfaces,  and  vortices  form  within  the  developing  plume  and  separation  occurs 
at  the  weapon  exit  plane  (see  Figure  1).  Further,  the  physical  properties  of 
the  plume  gas  are  often  sufficiently  different  from  those  of  air  that  it  is  of 
interest  to  take  this  into  account.  This  report  describes  a  model  for  such  a 
flow  based  upon  the  Ruler  equations  for  a  mixture  of  two  gases.  Marten's 
Total  Variation  Diminishing  (TVD)  scheme  (ref  1)  is  used  as  the  equation 
solver.  For  verification,  the  classic  problem  of  shock  diffraction  around  a 
corner  is  computed  and  compared  with  the  experimental  data  of  Skews  (refs 
2,3).  In  Reference  4,  the  model  is  used  to  analyze  the  origin  of  secondary 
shock  waves  frequently  observed  in  the  blast  signature  of  shoulder  fired 
weapons. 

THE  RULER  EQUATIONS 

Written  in  conservation  form,  the  Euler  equations  Lake  the  form 
iJU  d  F(U)  3G(U) 

—  +■ - 4- - 4-  W(U)  =  0  (1) 

at  ax  ay 


^Harten,  A.,  "High  Resolution  Schemes  for  Hyperbolic  Conservation  Laws," 
Journal  of  Computational  Physics,  Vol.  49,  1983,  pp.  357-393. 

-Skews,  8.  W. ,  "The  Shape  of  a  Diffracting  Shock  Wave,"  Journal  of  Fluid 
Mechanics ,  Vol.  29,  Part  2,  1967,  pp.  297-304. 

^Skews,  B.  W.,  "The  Perturbed  Region  Behind  a  Diffracting  Shock  Wave,"  Journal 
of  Fluid  Mechanics,  Vol.  29,  Part  4,  1967,  pp.  705-719. 

^Carofano,  G.  C.,  "Secondary  Waves  From  Nozzle  Blast,”  U.S.  ARDC  Technical 
Report  No.  ARLCB-TR-34028 ,  Benet  Weapons  Laboratory,  Watervliet,  NY,  October 
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S  »  S'/p0,  X  «  X'/R,  Y  *  Y'/R,  and  c  =  t' (P0/Po)^^/R  where  the  primed 
quantities  are  the  dimensional  counterparts  of  those  defined  above. 


The  pressure  is  related  to  the  thermodynamic  state  variables  p  and  e 
through  the  equation  of  state.  Rather  than  carrying  the  details  of  this 
relationship  through  the  algebra  of  determining  the  eigenvalues  and  eigen¬ 
vectors  of  the  Euler  equations,  the  equation  of  state  will  be  written  in  the 
general  form 

P  =  P(P,m,n,E,S)  =  P(p,e,Mf)  (3) 

where  Mf  Is  the  mass  fraction  of  species  "a”  defined  as 

Mf  =  S/p  (4) 


The  acoustic  speed,  c,  is  determined  from  the  equation  of  state  by 

2  .  3P1 
C  =  --J 

3p  s,Mf 

where  s  is  the  specific  entropy.  Using  the  following  general  relationship  for 
partial  derivatives  _ 

--)  f  --) 

P,Mf  3p  s,Mf  3P  e,Mf 

and  the  thermodynamic  identity  (ref  5) 


3p 


s,Mf 


3P 

3e 


3e 

P  =  ~) 

8p  s , Mf 


the  acoustic  speed  becomes 


P  3p, 


3P, 


c2  „  „)  f  -) 

p  3e  P,Mf  3p  e,Mf 


(5) 


50bert,  E.  F . ,  Concepts  of  Thermodynamics.  McGraw-Hill  Book  Co.  Inc.,  New 
York,  I960,  Equation  11-lb. 
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The  subscript  "Mf”  Is  used  to  Indicate  that  the  partial  derivatives  are 
evaluated  for  a  specific  mixture. 


In  what  follows,  other  partial  derivatives  of  pressure  appear  and 
useful  relationships  exist  among  them.  From  Eq.  (3), 

3P. 


some 


3p 


3P  dp  3e  dP  3Mf, 


ra 


,n,E,S  3P  e,Mf  3e  p,Mf  3P  ra,n,E  3Mf  p,e  3^  S 


3P  3P  3e, 

Pm  -  -")  =  ") 

dm  P,n,E,S  3e  p,Mf  dm  P,n,E,S 


3  P  3P ,  3e 

pn  -  ")  =  ~)  “) 

dn  P,m,E,S  3e  P,Mf  dn  P,ra,n,S 


3p 

pE  -  “) 


2)  2) 


3e  p,m,n,S  3e  P,Mf  3E  P,m,n,S 


3p  3p  3Mf 

Pg  = - j  = - ]  - ) 

3s  p,m,n,E  3%  p>e  3s  P 


Equations  (2)  through  (6)  can  be  combined  to  give 

pm  +  u  PE  *  0 
Pn  +  v  PE  =  0 


c2  *  Pp  +  MfPs  +  (U  -  u2  -  v2)Pe 
where  H  is  the  total  enthalpy  per  unit  mass,  defined  as 

E  +  P 


11 


(6a] 


(6b] 


(6c; 


(6d; 


(6e; 


(7: 

(8) 

(9) 


(10) 


The  matrix  manipulations  required  to  determine  the  eigenvalues  and 
eigenvectors  are  straightforward  and  will  not  be  presented  here.  Suffice  it 
to  say  that  the  combinations  of  variables  and  derivatives  in  Eqs.  (7)  through 
(10)  appear  frequently  and  lead  to  the  compact  and  recognizable  forms  given 


below* 


denote  the  Jacobian 

matrix 

3F(U)/3U, 

or 

0 

1 

0 

0 

0 

-u2  +  Pp 

2u  +  Pm 

Pn 

Pe 

Ps 

A  » 

-uv 

v 

u 

0 

0 

-uH  +  uPp 

H  +  uPm 

vPn 

u(l+PE) 

uPs 

-uMf 

Mf 

0 

0 

U 

(U) 


Setting  the  determinant  |A-axl|  equal  to  zero,  where 

matrix,  and  solving  for  the  vector  ax  of  eigenvalues 

,12  3  4  5 

(a  ,  a  ,  a  ,  a  ,  a  )  »  (u-c,  u,  u, 
X  X  X  X  X 


I  Is  the  Identity 

gives 

u+c,  u) 


(12) 


*The  author  began  this  study  based  on  a  report  written  In  early  1982  by  Harten 
at  the  Courant  Institute,  New  York  University.  That  report,  dealing  mostly 
with  theory,  eventually  appeared  as  Reference  l  In  1983.  Subsequently,  two 
more  reports  were  written  by  Yee,  Warming,  and  Harten  (refs  6,7)  which 
emphasized  applications  of  the  theory.  Thus,  as  an  aid  to  the  reader,  the 
present  report  will  conform  to  the  notation  of  the  later  works  where 
convenient,  but  some  of  the  symbols.  Indices,  and  even  the  "arrangement’*  of 
the  arrays  follow  the  preliminary  work  performed  by  the  author  In  developing 
his  computer  program. 

^Harten,  A.,  "High  Resolution  Schemes  for  Hyperbolic  Conservation  Laws," 
Journal  of  Computational  Physics,  Vol.  49,  1983,  pp.  357-393. 

^Yee,  H.  C.,  Warming,  R.  F.,  and  Harten,  A.,  "A  Hlgh-Resolutlon  Numerical 
Technique  For  Invlscld  Gas-Dynaralc  Problems  With  Weak  Shocks,"  Preprint  for 
Proceedings  of  the  Eighth  International  Conference  on  Numerical  Methods  in 
Fluid  Dynamics ,  Aachen,  West  Germany,  June  1982. 

?Yee,  H.  C.,  Warming,  R.  F.,  and  Harten,  A.,  "Implicit  Total  Variation 
Diminishing  (TVD)  Schemes  for  Steady-State  Calculations,”  NASA  Technical 
Memorandum  84342,  March  1983. 
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The  right-eigenvectors  Rx  ■  (R  •  R  »  R  »  R  »  R  )  are 

X  X  X  X  X 


1 

1 

0 

1 

u-c 

a 

0 

u+c 

V 

V 

1 

V 

H-uc 

H-c2/PE 

V 

11+uc 

Mf 

Mf 

0 

Mf 

1  2  3  4  5 

and  the  left-eigenvectors  =  (L  ,  L  ,  L  ,  L  ,  L  )  are 

X  X  X  X  X 


(u/c  +  r)/2 

1-r 

-v 

(-u/c  +  r)/2 
-Mf 


— (1/c  +  uq)/2  -vq/2 

uq  vq 

0  1 

(1/c  -  uq)/2  -vq/2 

0  0 


-rq/2 

rq 

0 

-rq/2 

i 


where 


r 

=*  'ps/pE 

(15) 

.  *  .  •  .  * 

*  '  •  .  ■  < 

q  =  Pg/c2 

»  r  = 

Pp/c2  =  l-q(H-u2-v2-Mf  f) 

(16) 
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the  Jacobian 

matrix 
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-uv 
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u 

0 
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-v2+Pp 

pra 

2vf?n 

PE 
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THE  ALGORITHM 


Operator  splitting  was  used  to  solve  Eq.  (I)  on  a  uniform  grid.  Harten's 

TVD  scheme  (refs  1,6,7)  was  applied  to  the  two  flux  terms  and  Euler's 

predictor-corrector  method  was  used  to  handle  the  source  term.  The  solution 
n 

at  time  t,  j,  was  advanced  to  time  t  +  2 At  using  the  following  sequence  of 
operations: 


n+2  n 

°l,j  *  LX  LY  Ls  lS  lY  ul,j  (21) 

where 

LX:  Ui, j  -  Ui,j  “  ”  (Fl+i/2,j  ”  Fi-i/2,j)  (22a) 

**  *  at  ** 

Ly:  Ulfj  =  Uifj  -  --  (G1>j+i/2  -  Gi^j-x/2)  (22b) 

—  **  ** 

Ui.j  =  Ulfj  -  At  mit5) 

LS: 

n+1  ** 

Ui.j  =  [Ui.j  +  u£>j  -  at  W(Uij)]/2 

a 

The  flux  Fi+l/2,j  Is  given  by 

A  5  ^  k 

Fl+l/2,j  *  [F(Ui,j)  +  F(U1+1>j)  +  (AX/At)  l  &i+i/2, j  Ri+l/2,j]/2 

k=l 

k  k  k  k  k  k 

® i+l / 2 , j  =  <81, j  +  81+1, j)  -  Q( vi+l/2 , j  +  ^i+l/2, j)  aL+l/2, j 


(22c) 

(23a) 

(23b) 


^-Harten,  A.,  "High  Resolution  Schemes  for  Hyperbolic  Conservation  Laws," 
Journal  of  Computational  Physics,  Vol.  49,  1983,  pp.  357-393. 

^Yee,  H.  C.,  Warming,  R.  F.,  and  Harten,  A.,  "A  High-Resolution  Numerical 
Technique  For  Invlscld  Gas-Dynamic  Problems  With  Weak  Shocks,”  Preprint  for 
Proceedings  of  the  Eighth  International  Conference  on  Numerical  Methods  In 
Fluid  Dynamics ,  Aachen,  West  Germany,  June  1982. 

^Yee,  H.  C.,  Warming,  R.  F.,  and  Harten,  A.,  "Implicit  Total  Variation 
Diminishing  (TVD)  Schemes  for  Steady-State  Calculations,"  NASA  Technical 
Memorandum  84342,  March  1983. 


(23c) 


Sl,j  *  3 1+1/2 , j  max[9,min( lat+l/2, jl »  al-l/2,j  8l+l/2,j]/8 
k  k 

sl+l/2, j  "  *l3*»(«1+l/2tj) 


k 

Yi+l/2,j 


k  k  k 

<Sl+l,j  -  81, j)/al+l/2, j 


k 

al+l/2 , j  *  0 


k 

al+l/2,j  =  0 


(23d) 


(23e) 


k  k 

v 1+1/2 , j  =  (Ax/iX)  at+i/2,j  ( 2 3 f ) 


Q(z)  *  z2  +  1/4 


(23g) 


k  k 

where  ai+i/2,j  and  Rt+l/2,j  are  given  In  Eqs.  (12)  and  (13),  respectively. 

The  subscript  1+1/2  signifies  that  all  of  the  values  In  the  functions  are  to 
be  evaluated  at  some  average  state  (Uj^j,  Ui+l^j).  In  this  study,  this  state 
was  taken  to  be 


u  =»  (ulf  j  +  uiKf  j)/2 

(24a) 

v  -  (v1>:J  +  v1+l>j)/2 

a 

(24b) 

c  =  (Cl>j  +  c1+l>j)/2 

(24c) 

11  ”  <Hl,j  +  lIl+l,j)/2 

(24d) 

Mf  =*  (Mf  lf  j  +  Mf  i+l ,  j)/2 

(24e) 

PE  “  (pEl ,  j  +  pCi-+l ,  j>/2 

(24f) 

r  “  (rl, j  +  ri+l, j)/2 

(24g) 

k 

The  vector  Is  computed  as  follows: 


*  pi+l,j  ~  pi,j  (25a) 

6ra  »  »>t+i ,  j  ~  mi ,  j  (25b) 

6n  *  ni+l,j  “  nl, j  (25c) 
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[• 


ft 

ft 

9 

9 

9 


-  r6p  +  q( 6G  -  ufim  -  v5n  -  I'6S) 


(26a) 


C2  =  (6m  -  u  6  p)  /  c 

(26b) 

ot1  =  (C,  -  C?)/2 

i+l/2,j 

(26c) 

a2  =6 p  -  Ci 

1+1/2 ,  j 

(26d) 

=  6n  -  v6p 

i+1/2, i 

(26e) 

a4  =  (Ci  +  C2)/2 

i+l/2,1 

( 2  6  f ) 

a-*  =  6g  -  Mf6p 

i+1/2,} 

(26g) 

The  flux  Gitj+i/2  ls  given  by  expressions  similar  Co  Eqs.  (23),  only  with 

k  k 

the  subscript  j  varying;  a^j+1/2  an<*  *ii,j+l/2  are  given  by  Eqs.  (18)  and 

k 

(19),  respectively.  The  vectoc  ^  j-4-1  /2  uses  analogues  of  Eqs.  (24)  and 
(25)  with  j  varying.  Finally 


Cl  - 


rop  +  q(5E  -  u6m  -  v6n  -  T6S) 

(27a) 

C2  =  ( 6n  -  v6p)/c 

(27b) 

a1  -  (Ci  -  C2 ) / 2 

i, j+1/2 

(27c) 

u2  =  3p  -  Ci 

i,j+l/2 

(27d) 

aJ  =  om  -  u6p 

i >  j+l/2 

(27e) 

u4  =  (Ci  +  C2)/2 

i, j+l/2 

(  2  7  f  ) 

a5  *  5S  -  ?lf  6p 

J  l  1  1  /  ^  1 

(27g) 

The  time  step,  At,  was  calculated  from  the  expression 


At 


0.85  AX 

maxi, jl®ax( lu^ j! , |vlf j|)  +  cj^j] 


all  i,j. 


The  CFL  factor  of  0.85  was  suggested  in  Reference  1. 

For  axisymmetric  calculations  (e  =  1),  the  coefficient  (n/Y)  in  the 
vector  W(u),  Eq.  (1),  becomes  indeterminant  on  the  axis  because  the  radial 
component  of  momentum  and  Y  both  vanish.  To  overcome  this  difficulty, 
L'Hospital's  rule  is  applied  to  Eq.  (1).  Since  all  of  the  terms  are  well 
behaved  at  Y  =  0  except  the  coefficient  (n/Y),  only  this  factor  is  changed  in 
the  limiting  process.  It  becomes 

3  n/  3Y  3n 

1  3Y 


The  following  second  order  accurate  Lagrange  differentiation  formula  is  used: 

3no  1 

- - - (-n_i  +  ni) 

3Y  2  AY 

where  the  n__^ ,  n0,  and  n^  are  the  three  points  at  j  «  -1 ,  0,  and  1, 
respectively.  By  symmetry  u_^  =  -ni ,  therefore, 

3n0  nx 
3Y  AY 


COMMENT:  The  reader  familiar  with  Marten's  scheme  will  recognize  that 
the  "artificial  compression"  terms  have  been  omitted  from  the  flux  expression 
in  Lq.  (23).  These  terms  were  found  to  be  quite  useful  for  "squaring  the 
corners"  at  shocks  and  contact  surfaces  but  frequently  caused  an  entropy 
violation  in  regions  where  a  shock  was  followed  by  an  expansion  or  where  an 
expansion  was  terminated  by  a  shock.  Reducing  the  amount  of  artificial 


*Harten,  A.,  "High  Resolution  Schemes  for  Hyperbolic  Conservation  Laws," 
Journal  of  Computational  Physics,  Vol.  49,  1983,  pp.  357-393. 


compression,  as  suggested  In  Reference  7,  would  alleviate  the  problem  In  one 
circumstance  only  to  have  it  reappear  In  another.  Because  such  regions  are 
common  to  the  problems  of  Interest,  the  decision  was  made,  quite  reluctantly, 
to  omit  the  compression  terms  over  the  entire  flow  field. 

Secondly,  the  reader  will  note  that  the  "average  state”  defined  above 
differs  from  the  two  choices  suggested  In  References  1,  6,  and  7.  One  choice 
Roe's  method,  uses  a  "mean  value  Jacobian”  such  that 

F(Ui+ifj)  -  F(Ui>j)  =  A(Ut>j,lJ1+1>j)<Ui+l  j  -  Ulfj) 

For  an  ideal  gas  It  Is  not  difficult  to  calculate  the  proper  u,  v, . . .  from 
this  expression,  but  for  more  general  equations  of  state  the  required  result 
Is  difficult  to  obtain. 

The  second  choice,  (Uj^j  +  Ul+llj)/2,  occasionally  produced  poor  results 
when  used  with  more  general  equations  of  state.  It  first  averages  the 

AAA 

conserved  variables,  then  calculates  the  required  quantities,  u,  v,  c,...  froi 
the  average  state. 

The  method  used  here,  Sq .  (24),  evaluates  the  variables  at  each  location 
then  forms  the  required  average.  This  worked  consistently  well  for  various 
state  equations;  for  an  ideal  gas,  It  produced  essentially  the  same  solutions 
as  were  obtained  with  Roe's  averaging. 

^Harten,  A.,  "High  Resolution  Schemes  for  Hyperbolic  Conservation  Laws,” 
Journal  of  Computational  Physics,  VoL.  49,  1983,  pp.  357-393. 

^Yee,  H.  C.,  Warming,  R.  F.,  and  Marten,  A.,  "A  High-Resolution  Numerical 
Technique  For  Invlscid  Gas-Dynamic  Problems  With  Weak  Shocks,”  Preprint  for 
Proceedings  of  the  Eighth  International  Conference  on  Numerical  Methods  in 
Fluid  Dynamics,  Aachen,  West  Germany,  June  1982. 

^Yee,  H.  C.,  Warming,  R.  F.,  and  Harten,  A.,  "Implicit  Total  Variation 
Diminishing  (TVD)  Schemes  for  Steady-State  Calculations,"  NASA  Technical 
Memorandum  84342,  March  1983. 


THE  STATE  EQUATION 


The  present  model  was  developed  as  an  aid  for  Interpreting  some 
experimental  results  obtained  with  a  blast  simulator  (ref  4).  Mixing  of  the 
helium  driver  gas  with  air  occurs  In  the  system  and  an  equation  of  state  for 
the  mixture  Is  needed. 

The  state  variables  In  the  Euler  equations  are  the  density,  p,  and  the 
specific  internal  energy,  e.  The  latter  has  units  of  energy  per  unit  mass, 
but  the  properties  of  a  mixture  of  gases  are  best  determined  on  a  molal  basis 
(ref  8).  Molal  units  will  therefore  be  used  in  the  Intermediate  steps  leading 
to  the  desired  function  P  ■  P(p,e,Mf). 

Consider  a  mixture  of  two  ideal  gases  containing  mass  of  species  "a** 

and  mass  of  species  "b".  The  mass  fraction,  Mf,  of  species  "a**  Is  then 

“a 

„f - 

ma  4-  n»b 


The  number  of  moles,  n,  of  each  species  Is 

raa  i»b 

na  *  .  %  “  ~ 

Ma  % 

where  Ma,  Mb,  are  the  respective  molecular  weights.  The  mole  fractions  are 


"b 

,  xb  - - - - 


For  the  mixture 


na  +  nb  na  +  nb 

M  -  xaMa  ♦  xbMb 


^Carofano,  G.  C.,  "Secondary  Waves  From  Nozzle  Blast,"  U.S.  ARDC  Technical 
Report  No.  ARLCB-TR-84028 ,  Benet  Weapons  Laboratory,  Watervllet,  NY,  October 
1984. 

®Obert,  E.  F.,  Concepts  of  Thermodynamics.  McGraw-Hill  Book  Co.  Inc.,  New 
York,  1960,  Chapter  8 .  ” 


CV 


''a'-va  xbcvb 


cp  =  xacpa  xbcpb 

where  cv  and  Cp  are  the  molal  specific  heats  at  constant  volume  and  constant 
pressure,  respectively.  For  each  gas 


-va 


*o  *o 

Z-i  •  c'rb  ■  7h:r 

where  Rq  Is  the  universal  gas  constant  In  molal  units.  The  specific  Internal 
energy  of  the  mixture  is 

(xacva  ■*"  xbcvb)T 


where  T  Is  the  temperature;  the  reference  state  Is  taken  as  absolute  zero. 

The  state  equation  for  the  mixture  Is 

-  PR0T 

P  - - 

M 

or,  using  the  above  definitions, 

Pe 

P  - -  (28) 

xa  xb 

( - «• - ) 

<b~l 

The  specific  heat  ratio  for  the  mixture  Is  defined  as 

y  -  Cp/ cv 

With  a  little  algebra  the  following  results  are  obtained: 


1  xa  xb 

Y-l  Ya-1  Yb-1 

MfOa  "  °fb>  +  °fb 

Y  - - 

Mf ( l-o)  +  o 

3  Y  °(  fa  -  *b> 

Mf  3Mf  [Mf(l-O)  +  o]2 


(29) 


(30) 


(31) 
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Ma(Ya-i> 

0  - — - 

MbC^b"1) 


(32) 


P  *  (  Y-l )  Pe 

(33) 

PE  -  (Y-l) 

(34) 

ps  ■  e,Mf 

(35) 

c2  =  YP/p 

(36) 

-eY 

Mf 

r  „  - 

(37) 

(Y-l) 

where  Eqs.  (5),  (6),  and  (15)  have  been  used  to  obtain  the  last  four 
equations.  The  specific  Internal  energy  is  calculated  from  the  conserved 
variables  using  Eq.  (2). 

For  a  given  pair  of  gases  Ya,  Yb,  Ma  and  Mb  are  constants  along  with  a  In 
Eq.  (32).  Therefore,  only  Eqs.  (30),  (31),  and  (33)  through  (37)  are  actually 
used  In  the  algorithm  given  in  the  previous  section. 

When  Ya  3  Yb  and  Ma  »  Mb  all  of  the  results  reduce  to  those  given  in 
References  1,  6,  and  7,  and  the  last  of  Eqs.  (1)  for  S  is  not  needed. 

4owever,  as  will  be  shown  below,  even  for  a  single  gas  species  some  useful 
Information  can  often  be  obtained  by  retaining  it. 


^Harten,  A.,  "High  Resolution  Schemes  for  Hyperbolic  Conservation  Laws," 
Journal  of  Computational  Physics,  Vol.  49,  1983,  pp.  357-393. 

^Yee,  H.  C.,  Warming,  R.  F.,  and  Harten,  A.,  "A  Htgh-Resolutlon  Numerical 
Technique  For  Inviscid  Gas-Dynamic  Problems  With  Weak  Shocks,”  Preprint  for 
Proceedings  of  the  Eighth  International  Conference  on  Numerical  Methods  In 
Fluid  Dynamics,  Aachen,  West  Germany,  June  1982. 

^Yee,  H.  C.,  Warming,  R.  F.,  and  Harten,  A.,  "Implicit  Total  Variation 
Diminishing  (TVD)  Schemes  for  Steady-State  Calculations,”  NASA  Technical 
Memorandum  34342,  March  1983. 


BOUNDARY  CONDITIONS 

Along  the  X-axls  the  symmetry  condition  was  applied  in  the  usual  manner 
by  using  mirror  images  of  points  in  the  active  mesh  but  with  the  sign  of  the 
v-component  of  velocity  changed.  The  Inflow  boundary  was  handled  either  by 
using  a  known  solution,  as  in  the  test  problem  of  the  next  section,  or  by 
zeroth-order  extrapolation  of  characteristic  variables  (refs  9,10)  as  will  be 
explained  below. 

For  the  X-sweep  along  the  vertical  wall  (see  Figure  l),  the  reflection 
method  was  used.  For  the  Y-sweep  flow  separation  was  assumed  and  an  extra 
point  was  added  at  the  corner  -  point  “b"  in  the  sketch  below.  The  wall 
boundary  condition  dictates  that  u  -  0.  Since  there  can  be  no  flow  across  the 
contact  surface,  v  =  0  must  apply.  Therefore,  the  same  particle  remains  there 


separated 

boundary 


flow 


^Uarten,  A.,  "On  a  Large  Time-Step  High  Resolution  Scheme,"  ICASE  Report  No. 
82-34,  NASA  Langley  Research  Center,  Hampton,  VA,  November  15,  1982. 
l^Yee,  H.  C.,  Beam,  R.  M.,  and  Warming,  R.  F.f  "Boundary  Approximations  for 
Implicit  Schemes  for  One-Dimensional  Inviscld  Equations  of  Gas  Dynamics," 
AIAA  .Journal,  Vol.  20,  No.  9,  September  1982,  pp.  1203-1211. 
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and  Its  mass  fraction  remains  constant,  l.e.,  Mf  ■  1  if  the  particles  In  the 
environment  are  Identified  as  species  "a".  The  remaining  unknowns  are  density 
and  energy  so  two  equations  are  needed.  Because  v  ■  0  at  the  corner, 
information  should  reach  It  from  point  "a"  only  through  the  (v-c)  character¬ 
istic  variable.  The  characteristic  variables  for  the  Y-sweep  are  w  *  LyU 
where  by  is  given  by  Eq.  (20).  The  first  characteristic  variable  correspond¬ 
ing  to  the  (v-c)  eigenvalue  is 

AAA  .X  A  A  A  A  A  A  A 

w1  =  [(v/cl-r)p  -  uqm  -  (1/c  +  vq)n  +  qE  -  TqS)/2 

The  approximation  w1  =  w1  is  used  at  time  t  +  At  with  P,  c,  ...  evaluated  at 
b  a 

the  average  of  the  two  states. 

The  second  characteristic  variable  corresponding  to  the  (v)  eigenvalue 
is  assumed  to  remain  constant.  It  is  given  hy 

A  AA  A  A  A  AA 

w2  “  (l-r)p  +  uqm  +  vqn  -  qE  +  TqS 

The  assumption  Is  that  w2  at  time  t  is  equal  to  w2  at  t  +  At.  The  average 

..  b  a  b 

state  values  of  u,  v,  and  f  are  always  zero. 

The  procedure  is  as  follows: 

1.  Apply  the  operator  by  to  the  points  along  the  wall  starting  with 

point  "a”.  This  requires  Information  at  point  “b”  at  time  T,  which  is 

available,  and  Information  at  a  second  point  below  "b"  because  Marten* s  method 

k 

is  a  five-point  scheme;  it  is  needed  to  calculate  at  "b" .  Since  it  is 

k 

not  available,  8i,j  at  point  "a"  is  used. 

2.  Solve  the  two  characteristic  relations  listed  above  to  obtain  Pb  and 
Eb  at  time  x  +  At.  Iteration  is  required  because  the  average  state  values 
depend  on  the  new  solution  at  ”b”. 


One  final  note:  It  was  found  through  experience  that  the  X-sweep 
generally  could  not  be  applied  to  the  first  one  or  two  points  along  the  wall 
without  generating  physically  unrealistic  results  there.  It  will  be  seen  in 
the  next  section  that  this  extremely  crude  handling  of  the  corner  flow  does 
not  seem  to  have  a  drastic  Impact  on  the  rest  of  the  flow  field.  This  may  be 
true  because  the  inflow  boundary  condition  is  so  dominant. 

A  TEST  PROBLEM 

The  classic  shock  diffraction  experiments  of  Skews  (refs  2,3)  contain 
many  of  the  features  of  the  flow  sketched  in  Figure  1  and  afford  an  excellent 
test  problem.  Superimposed  upon  the  density  contour  plot  in  Figure  2  are 
Skews'  data  for  a  planar  shock  (e  =  0)  with  Mach  number  Mq  *  3  diffracting 
around  a  90  degree  corner.  The  flow  behind  the  shock  is  supersonic  for  this 
case  so  the  disturbed  region  lies  entirely  downstream  of  the  corner.  The 
inflow  boundary  conditions  upstream  were  constant  during  the  calculation. 

The  construction  of  Figure  2  is  facilitated  by  the  self-similar  nature  of 
the  planar  flow  (ref  11).  Thus,  the  unsteady  numerical  solution  can  be 
compared  to  the  unsteady  laboratory  experiment  at  arbitrary  times  if  the  data 
are  scaled  using  the  similarity  variables 

x'  X  y'  Y 

co«-  /yx  cQt  /yx 


^Skews,  3.  W. ,  "The  Shape  of  a  Diffracting  Shock  Wave,”  Journal  of  Fluid 
Mechanics,  Vol.  29,  Part  2,  1967,  pp.  297-304.  -  - 

^Skews,  8.  W. ,  "The  Perturbed  Region  Behind  a  Diffracting  Shock  Wave,” 
Journal  of  Fluid  Mechanics,  Vol.  29,  Part  4,  1967,  pp.  705-719. 

Hjones,  P.  M. ,  Martin,  P.  M.  E.,  and  Thornhill,  C.  K. ,  "A  Note  on  the 

Pseudo-Stationary  Flow  Behind  a  Strong  Shock  Diffracted  or  Reflected  at  a 
Corner,"  Proceedings  of  the  Royal  Society,  A209 ,  1959,  pp.  238-248. 


Here  cn  Is  Che  acoustic  speed  In  the  undisturbed  environment;  for  an  Ideal  gas 
co  *  / y’p07"p0.  The  factor  /y  appears  In  the  denominator  because  t  Is  non- 
dlmenslonallzed  with  *^7 PQ  In  the  present  study.  The  numerical  data  In 
Figure  2  were  obtained  after  135  cycles  on  a  uniform  grid  with  150  cells  In 
the  X-dlrectlon  and  190  cells  In  the  Y-dlrection;  Y  was  taken  as  1.4  for  air. 

The  predicted  shock  position  agrees  with  Skews'  measurements  (see  Table 
l)  to  within  three  percent  everywhere.  Note  that  the  Inflection  point  In  the 
experimental  data  near  the  wall  Is  also  present  In  the  numerical  solution.  It 
corresponds  to  a  Mach  reflection. 

TABLE  l.  SHOCK  POSITION  (ESTIMATED  FROM  FIGURE  3  OF  REF  2) 


£ 

n 

2.75 

0.00 

2.41 

0.65 

1.99 

1.15 

1.48 

1.48 

0.95 

1.65 

0.44 

1.64 

0.00 

1.58 

The  disturbance  which  propagates  downward  into  the  uniform  flow  behind 
the  shock  consists  of  two  segments,  the  Mach  line  or  leading  edge  of  the 
Prandtl-Meyer  expansion  centered  at  the  corner  and  a  circular  soundwave. 
Using  expressions  given  by  Skews  (ref  2)  and  the  nomenclature  of  Figure  3 


^Skews,  B.  W.,  "The  Shape  of  a  Diffracting  Shock  Wave,”  Journal  of  Fluid 
Mechanics ,  Vol.  29,  Part  2,  1967,  pp,  297-304. 


the  following  can  be  written: 


cit  ci 


2(  Y-l) 
( Y+l ) 2 


2  2 
(Mq-IKY+I/Mq) 


- (M0 

(Y4-1) 


«1  Ul  Cq  Up 

Ml  -  ~  - - -  “ 

C1  c0  Cl  < 

M  =  sin" ^(1/Mi) 

2  3 

tan2  a  =  Up[(Y-l)Mo/2  +  U/Mo 


^O  ~ 

n  -  -C0 

tan  a 

Cl  -  Up 

hi  *  -te 

Ct  *  Ci  -  <  sin  p 

\  *  -K 

cos  U 

where  <  is  the  ratio  of  acoustic  speeds  across  the  shock  and  up  Is  the 
dimensionless  particle  velocity  behind  the  shock.  Points  along  the  Mach  line 
are  computed  from 

n  =  -C  tan  u 

and  intermediate  points  on  the  circular  arc  are  given  by 

(C-Cx)2  +  n2  =  n2 

The  numerical  solution  smears  the  leading  edge  of  the  Prandtl-Meyer 
expansion  slightly  so  that  the  first  contour  line  does  not  exactly  coincide 
with  the  Mach  line.  The  corner  approximation  may  also  have  influenced  the 
result.  Elsewhere  the  agreement  is  satisfactory. 

The  terminator  angle,  6,  was  estimated  to  be  20  degrees  ((ref  3),  Figure 


^Skews,  B.  W. ,  "The  Perturbed  Region  Behind  a  Diffracting  Shock  Wave,"  Journal 
of  Fluid  Mechanics,  Vol.  29,  Part  4,  1967,  pp.  705-719. 


6).  It  marks  the  end  of  the  expansion  fan  and  the  beginning  of  a  small  region 
of  uniform  flow  as  discussed  in  Reference  11.  The  contours  in  Figure  2 
conform  to  this  description. 

The  slipstream  angle,  u>,  was  estimated  to  be  40  degrees  ((ref  3),  Figure 
5).  Skews  used  it  as  a  reference  line  along  which  was  measured  the 
intersections  of  the  second  (recompression)  shock  and  the  contact  surface. 

The  contact  surface  in  Figure  2,  indicated  by  the  thin  solid  line  following 
the  shock,  represents  the  boundary  separating  the  particles  which  were 
processed  by  the  shock  before  it  reached  the  corner  from  those  processed  after 
it  diffracted.  it  is  obtained  from  the  numerical  solution  by  “tagging'’  one 
set  of  particles  as  "species  b”  and  assigning  them  a  mass  function  Mf  =  0. 

Both  sets  of  particles  are  given  the  same  molecular  weight  and  specific  heat 
ratio  so  T  is  zero  everywhere.  The  solution  evolves  independently  of  the 
species  equation.  The  contact  surface,  initially  a  vertical  line  at  C  *  0,  is 
simply  convected  downstream  by  the  velocity  field.  Note  that  it  passes 
through  the  low  point  of  the  soundwave.  This  is  correct  because  all  of 

the  particles  below  the  soundwave  move  at  the  particle  velocity  behind  the 
shock,  Up  = 

The  intersection  of  the  contact  surface  with  the  slipstream  angle, 
indicated  by  the  asterisk  in  Figure  2,  was  estimated  to  be  1.65  ((ref  3), 
Figure  8c).  Its  coordinates  5C,  nc  are 

^Skews,  B.  W.,  "The  Perturbed  Region  Behind  a  Diffracting  Shock  Wave,"  Journal 
of  Fluid  Mechanics,  Vol.  29,  Part  4,  1967,  pp.  705-719. 

I  * Jones,  P.  M.,  Martin,  P.  M.  E.,  and  Thornhill,  C.  K. ,  "A  Note  on  the 
Pseudo-Stationary  Flow  Behind  a  Strong  Shock  Diffracted  or  Reflected  at  a 
Corner,"  Proceedings  of  the  Royal  Society,  A209,  1959,  pp.  238-248. 


Sc  =  1.65  cos  40  -  1.264 
nc  =  1.65  sin  40°  =■  1.061 

The  numerical  solution  comes  very  close  to  passing  through  this  point. 

It  should  be  noted  that  the  contact  surface  described  above  and  that 
discernible  in  a  shadowgraph  are  not  the  same.  The  contact  surface  discussed 
by  Skews  actually  represents  the  boundary  separating  those  particles  which 
were  processed  by  the  planar  shock  from  those  processed  by  the  diffracted 
shock.  It  is  associated  with  the  entropy  gradient  along  the  curved  portion  of 
the  shock  and  can  be  identified  In  Figure  4  as  the  entropy  contour  which 
intersects  the  point  where  the  shock  ceases  to  be  planar.  This  contour 
coincides  wiLh  the  "kink"  in  each  of  the  density  contours  in  Figure  2  which 
renders  it  visible  on  a  shadowgraph.  At  the  point  where  Skews'  measurement 
was  taken  the  entropy  contour  and  the  contact  surface  nearly  conicide  which 
accounts  for  the  agreement  noted  above. 

The  second  shock  location,  indicated  by  the  star  in  Figure  3,  was 
estimated  to  be  1.2  ((ref  3),  Figure  7c).  Its  coordinates  £s,  ns  are 

Cs  *  1.2  cos  40°  =  0.919 
ng  =  1.2  sin  40°  =  0.717 
Again,  the  numerical  solutton  does  fairly  well. 

The  vortex  angle,  <|»,  and  its  location  Indicated  by  the  diamond  in  Figure 
2,  were  estimated  to  be  48  degrees  and  1.15,  respectively  ((ref  3),  Figure  9). 
Its  coordinates,  5V,  nv  are 


^Skews,  B.  W. ,  "The  Perturbed  Region  Behind  a  Diffracting  Shock  Wave,"  Journal 
of  Fluid  Mechanics,  Vol.  29,  Part  4,  1967,  pp.  705-719. 


0.669 


£v  -  1.15  cos  48°  - 

-  1.15  sin  48°  -  0.743 

Judging  by  the  density  contours,  the  proper  vortex  angle  seems  to  have  been 
computed  but  its  position  appears  somewhat  closer  to  the  corner  than  measured 
by  Skews.  However,  Identifying  the  vortex  center  from  a  shadowgraph  involves 
a  greater  degree  of  uncertainty  than  the  other  measurements  so  the  apparent 
discrepancy  may  not  be  serious. 

CONCLUSION 

The  test  problem  of  the  last  section  demonstrates  that  Harten* s  TVD 
method  can  be  used  to  obtain  reasonably  accurate  solutions  for  the  flow 
problems  of  interest.  In  a  companion  report  (ref  4),  the  method  Is  applied  to 
a  problem  Involving  two  gases. 


^Carofano,  G.  C.,  "Secondary  Waves  From  Nozzle  Blast,"  U.S.  ARDC  Technical 
Report  No.  ARLCB-TR-84028 ,  Benet  Weapons  Laboratory,  Watervliet,  NY,  October 
1984. 
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